Conformal Prediction for Smartphone Telemetry

Group 46 - Antonio Pagnotta, Beatrice Mazzocchi, Domenico Mangieri

Abstract

Our work focus on leveraging Functional Data Analysis (FDA) and Conformal Prediction to objectively classify human physical activity,specifically Sedentary, Walking, and Running states,using asynchronous smartphone sensor telemetry. To transition from raw, noisy multi sensor data to a reliable analytical foundation, the pipeline first synchronizes all readings onto a universal 10Hz grid and isolates the rotation invariant linear acceleration magnitude. Rather than relying on rigid, preassigned experimental labels that might ignore human micro variations, the subject’s true physical state is dynamically segmented using B-spline smoothing and K-Means clustering applied to a 5 second rolling variance.

With the data cleanly partitioned and filtered for purity, the analysis shifts to an FDA approach, treating 5 second epochs as continuous mathematical curves to extract distinct biomechanical signatures. By evaluating functional means, variance, and the first derivative of acceleration (Jerk), we establish a clear structural understanding of each movement type. Building upon these functional representations, we subsequently apply a Conformal Prediction framework utilizing Functional Principal Component Analysis (FPCA) and Mahalanobis distances to establish, distribution free boundaries for normal behavior. By constructing 90% confidence ellipsoids in the latent space and translating them into time domain prediction bands, the model successfully hits its empirical coverage targets across all activities. Ultimately, this geometric evaluation reveals that low energy activities (like resting) generate tightly constrained regions capable of perfectly isolating micro movements and rejecting over 93% of anomalies, whereas higher energy activities naturally form larger, more encompassing ellipsoids due to their intrinsic biomechanical variance.

Data Preprocessing

Once obtained the data, first thing we do is transform the raw, asynchronously recorded sensor data, into a clean, synchronized, and evenly spaced dataset. Our dataset contains high-dimensional telemetry from multiple sources, including Compass, Ambient Light, and Magnetic Rotation sensors. By a quick look at the data, we noticed that the sensors rarely recorded their read value at the exact same millisecond, so this pipeline forces all sensor readings, regardless of their type, onto a shared timeline so they can be accurately compared and analyzed row by row.

In order to do this, we first order the raw data by timestamp and converts absolute time (milliseconds) into a relative, easy to read scale starting at exactly 0, establishing a universal 10Hz sampling rate and creating a target timeline with exactly one row every 0.1 seconds. While we keep the full sensor suite synchronized for a complete record, we primarily target the Linear Accelerometer data for the core of our analysis. We used the zoo package to calculate missing values. In this way, if a sensor didn’t record a value at an exact 0.1s interval, the script mathematically estimates it by drawing a straight line between the nearest surrounding data points.

df_raw <- read.csv("dati_sgravati (1).csv")

# Sort by timestamp and convert to relative seconds
df_raw <- df_raw[order(df_raw$timestamp), ]
t0 <- min(df_raw$timestamp)
df_raw$time_sec <- (df_raw$timestamp - t0) / 1000

sensor_cols <- setdiff(names(df_raw), c("timestamp", "time_sec"))

# Interpolate to common 10Hz grid
dt <- 0.1
t_grid <- seq(from = 0, to = max(df_raw$time_sec), by = dt)
df_aligned <- data.frame(time = t_grid)

interpolate_sensor <- function(time_col, val_col, new_grid) {
  valid <- !is.na(val_col)
  if (sum(valid) < 2) return(rep(0, length(new_grid)))
  z <- zoo(val_col[valid], time_col[valid])
  as.numeric(na.approx(z, xout = new_grid, rule = 2))
}

for (s in sensor_cols) {
  df_aligned[[s]] <- interpolate_sensor(df_raw$time_sec, df_raw[[s]], t_grid)
}

B-Spline Smoothing and Segmentation

In our experimental setup, the subject was instructed to first run, then remain sedentary, and finally walk. However, real human movement is rarely perfectly constant. During the “Running” interval, the subject might have briefly stopped to catch their breath (sedentary) or slowed down (walking). If we simply hard coded the entire first time block as “Running,” we would inadvertently label standing still as a “run,” corrupting all our subsequent functional data analysis. To solve this, we built a pipeline to objectively extract the true physical state at any given moment, rather than relying on what the subject was supposed to be doing.

To initiate this process, we first had to standardize the raw sensor data. Although our hardware captures rotation and orientation (via the Magnetic Rotation and Compass sensors), we have elected to use Linear Acceleration Magnitude as our primary analytical feature to ensure rotation invariance.

Next, to filter out the high frequency hardware jitter, we fit a continuous mathematical B-spline curve through the data. As a practical measure to prevent the system from crashing, this smoothing was applied iteratively in chunks. Once we had this clean signal, we calculated how much the movement fluctuated over a rolling 5 second window, operating on the premise that high variance indicates intense movement, while near zero variance indicates stillness.

This variance metric is exactly where we account for the micro variations in the user’s actual behavior. By applying K-Means clustering to our rolling variance, we automatically found the three true intensity clusters hidden inside the data. This means if the subject stopped for 10 seconds during the planned “Running” phase, K-Means organically detected the drop in variance and correctly reassigned those specific seconds to “Sedentary.”

Crucially, because this K-Means pre segmentation successfully divides our data into distinct, homogeneous classes, the resulting distribution of each class’s features is intrinsically unimodal. This strategic step makes the subsequent application of a Gaussian Mixture Model computationally redundant, opting for a Single Multivariate Gaussian.

Finally, as a last quality control step, we applied Label Filtering. We ran a median filter over the K-Means output labels to iron out physical impossibilities. For example, if the model classified a single 0.1 second frame as “Running” sandwiched between continuous seconds of “Sedentary,” the filter smoothed it out, ensuring our final activity blocks made chronological and physical sense.

# Smooth Acceleration magnitude to remove high-frequency jitter
if("LinearAccelerometerSensor" %in% names(df_aligned)) {
  raw_mag <- abs(df_aligned$LinearAccelerometerSensor)
} else {
  raw_mag <- sqrt(df_aligned$AccX^2 + df_aligned$AccY^2 + df_aligned$AccZ^2)
}

# Smoothing in chunks to save memory
chunk_size <- 2000
num_chunks <- ceiling(length(t_grid) / chunk_size)
smooth_mag <- numeric(length(t_grid))

for(i in 1:num_chunks) {
  idx <- ((i-1)*chunk_size + 1):min(i*chunk_size, length(t_grid))
  nbasis <- max(4, round(length(idx) / 4))
  basis_obj <- create.bspline.basis(rangeval = range(t_grid[idx]),
                                    nbasis = nbasis, norder = 4)
  fd_par <- fdPar(basis_obj, Lfdobj = 2, lambda = 0.01)
  sm_res <- smooth.basis(t_grid[idx], raw_mag[idx], fd_par)
  smooth_mag[idx] <- as.numeric(eval.fd(t_grid[idx], sm_res$fd))
}
df_aligned$SmoothMag <- smooth_mag

# Segment via Rolling Variance + K-Means
window_size <- round(5 / dt) # 5 second rolling window
roll_var <- rollapply(df_aligned$SmoothMag, width = window_size,
                      FUN = var, fill = 0, align = "right")

set.seed(123)
log_var <- log(roll_var + 1e-6)
km <- kmeans(log_var, centers = 3, nstart = 20)
centers_sorted <- sort(km$centers)

thresh1 <- exp(mean(centers_sorted[1:2]))
thresh2 <- exp(mean(centers_sorted[2:3]))

df_aligned$Activity <- cut(roll_var, breaks = c(-Inf, thresh1, thresh2, Inf),
                           labels = c("Sedentary", "Walking", "Running"))

# Clean labels with median filter to remove micro glitches
act_smooth <- runmed(as.integer(df_aligned$Activity), k = round(1/dt)+1)
df_aligned$Activity <- factor(act_smooth, levels=1:3,
                              labels=c("Sedentary", "Walking", "Running"))

FDA Exploration

Here we shift from standard row by row data manipulation into Functional Data Analysis (FDA). Instead of analyzing thousands of individual data points, we treat the chunks of time (5 second windows) as a single, continuous mathematical curve and by analyzing the shape of these curves, extract the true physical “signatures” of different activities.

To prepare the data for this analysis, we first apply a filtering step to chop the continuous timeline into manageable 5 second epochs. This is done to filter out “messy” windows (e.g., when the user is actively transitioning from standing to walking), keeping only the windows that are at least 80% pure in a single activity. Once these pure epochs are isolated, we proceed to FDA Object Creation, which converts the remaining discrete matrices back into continuous B-spline curves.

With our mathematical functions established, we can begin our visual exploration. We start by looking at every single 5 second curve on top of each other (spaghetti plots) . By forcing a shared Y-axis across the panels, we can instantly and visually compare the physical intensity of Sedentary vs. Walking vs. Running. To gain deeper biomechanical insights beyond raw magnitude, we performed a Jerk analysis, computing the first derivative of the acceleration. In physics, the rate of change of acceleration is called “Jerk,” and isolating this metric perfectly highlights sharp, sudden impacts, such as a foot hitting the pavement.

Shifting our focus from individual strides to overall trends, we check the functional means and variance as a statistical summary . We plot the “average” movement signature for each activity, wrapped in a shaded ribbon representing the standard deviation (to understand how much individual strides vary from the average).

Finally, to rigorously identify any anomalies within those general trends, we check the functional boxplots . This allows us to check the functional mean (most typical movement curve for that activity) and distinctly highlights any abnormal, irregular strides as outliers in black.

# Slice into 5 second fixed functional observations
pts_per_window <- round(5 / dt)
window_duration <- 5
n_windows <- floor(nrow(df_aligned) / pts_per_window)
trunc_len <- n_windows * pts_per_window

acc_matrix <- matrix(df_aligned$SmoothMag[1:trunc_len],
                     nrow = pts_per_window, ncol = n_windows)
activity_matrix <- matrix(as.character(df_aligned$Activity[1:trunc_len]),
                          nrow=pts_per_window)

get_majority_label <- function(v) {
  tab <- table(v)
  lab <- names(tab)[which.max(tab)]
  prop <- max(tab) / sum(tab)
  c(label = lab, prop = prop)
}

# Filter for pure activity windows
maj_stats <- apply(activity_matrix, 2, get_majority_label)
maj_labels <- maj_stats["label", ]
maj_props  <- as.numeric(maj_stats["prop", ])

keep <- maj_props >= 0.8  # Keep windows that are at least 80% one activity
acc_matrix <- acc_matrix[, keep, drop = FALSE]
window_labels <- factor(maj_labels[keep],
                        levels = c("Sedentary","Walking","Running"))

# Create FDA Object
window_time <- seq(0, window_duration, length.out = pts_per_window)
fda_basis <- create.bspline.basis(rangeval = c(0, window_duration),
                                  nbasis = 15, norder = 4)
acc_fd <- Data2fd(window_time, acc_matrix, fda_basis)


# Find global y axis limits so all 3 plots share the same scale
plot_grid_eval <- seq(0, window_duration, length.out = 100)
ylim_acc <- range(eval.fd(plot_grid_eval, acc_fd))

In these 5‑second FDA windows, Sedentary curves is flat near zero, confirming that the segmentation plus the purity filter isolates motionless periods. Walking windows show medium level, relatively regular oscillations, representing periodic steps with moderate variability. Running windows populate the entire plot with dense and variable trajectories, reflecting a more heterogeneous impact. Using a common y axis across panels makes the intensity ordering clear to visualize.

# 1st derivative
acc_deriv_fd <- deriv.fd(acc_fd, 1)

# Find global y axis limits for the derivative
ylim_jerk <- range(eval.fd(plot_grid_eval, acc_deriv_fd))

This jerk plot shows how sharply acceleration changes within each 5 second window. In Sedentary, curves are essentially flat around zero, confirming that the smoothing plus segmentation removed almost all irregular motion. Walking shows small, localized bursts at the end of the window, after a steady smooth motion. Running displays large, frequent oscillations over all the interval, indicating sustained high frequency impacts and much more irregular dynamics.

plot_grid <- seq(0, window_duration, length.out = 100)
eval_acc <- eval.fd(plot_grid, acc_fd)

# Build data frame with proper alignment
n_curves <- ncol(eval_acc)
df_eval <- data.frame(
  Time = rep(plot_grid, n_curves),
  Acc = as.vector(eval_acc),
  CurveID = rep(1:n_curves, each = length(plot_grid))
)

# Assign activity labels by matching curve indices to window_labels
df_eval$Activity <- rep(as.character(window_labels), each = length(plot_grid))

df_summary <- df_eval %>%
  group_by(Time, Activity) %>%
  summarise(
    Mean_Acc = mean(Acc, na.rm = TRUE),
    SD_Acc = sd(Acc, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(Activity = factor(Activity,
                           levels = c("Sedentary", "Walking", "Running")))

This plot summarizes the average acceleration profile of each activity and its variability over the 5 second window. Sedentary stays essentially flat near zero confirming almost no movement and very low variance. Walking has the highest mean level and widest band, reflecting greater variability despite moderate intensity. Running is slightly lower in mean and variance indicating more consistently high amplitude. The overlapping means tells us that class specific FPCA ellipsoids are needed to a deeper analysis.

These functional boxplots summarize the distribution of all 5 second curves using depth, so the central band represents the most typical trajectory for each activity. Sedentary has an extremely tight box and almost no black outlier curves, confirming very homogeneous idle behaviour after segmentation and purity filtering. Walking and especially Running show thicker boxes and a few clear outliers, corresponding to atypical strides or short transitions that deviate from the dominant pattern.

Conformal Prediction

Now we transition from descriptive statistics into the actual Conformal Prediction framework to quantify uncertainty. Our primary objective is not simply to explore which movements are similar, but to establish a rigorous statistical rule for defining “normal” activity signatures. For this reason, we prioritized the approach described in section 3 of the main reference paper over the purely exploratory clustering (pseudo densities and clustering trees) found in section 4.

The method in section 3 allows for the construction of simultaneous prediction bands, providing a strict statistical guarantee (controlled by conf_alpha) about how likely a classification is correct. These bands effectively create a “channel” or “tube” around our linear accelerometer signal. If a recording stays within this tube, it is confirmed as that specific activity with a confidence guarantee. However, if the signal falls outside these boundaries, the system can immediately identify it as an anomaly or a transition to a different activity, such as Shopping or Sedentary.

To construct these robust profiles (or ellipsoidal confidence regions), the first operation is Functional PCA (FPCA). By projecting our high frequency data into a lower dimensional space, we filter out sensor noise and focus on the structural movement patterns. This technique reduces complex 5 second curves into a few “Principal Components” (PCs). Instead of dealing with 50 distinct data points, a curve is now efficiently represented by just 2 or 3 numbers that successfully capture 90% of its overall shape.

Once the curves are reduced to these core components, we apply the Mahalanobis Distance. Unlike a standard Euclidean distance that treats all directions equally, Mahalanobis accounts for how the data is stretched, scaled, or correlated in our multi dimensional PC space. We use it to calculate a Conformity Score, essentially measuring how “weird” or “normal” a curve is relative to its assigned class. This is particularly effective because it correctly shapes our confidence regions into ellipsoids, naturally fitting the actual variance of the subject’s movements.

Using these conformity scores, we find a cutoff value such that our target confidence level (e.g. 90%) of our training data falls securely inside the boundary.

With this threshold established, we generate our prediction bands. This process converts the abstract PC thresholds back into the physical time domain, visualizing the aforementioned “envelope” or “tube.” If a future movement curve stays completely within this tube, it is statistically guaranteed to be consistent with that specific activity.

Finally, to ensure our theoretical guarantees hold up in practice, we check if 90% of our actual “Walking” data falls within the constructed “Walking” band, and verify if “Running” data correctly falls outside that very same “Walking” band.

activities <- c("Sedentary", "Walking", "Running")
conf_alpha <- 0.10

# Initialize data frames to collect plotting data
df_bands_all <- data.frame()
df_proj_all <- data.frame()
grid_eval <- seq(0, window_duration, length.out = 100)

# Diagnostic print to check labels
cat("Current distribution of functional windows:\n")
## Current distribution of functional windows:
print(table(window_labels, useNA = "always"))
## window_labels
## Sedentary   Walking   Running      <NA> 
##       405       211       448         0
coverage_results <- data.frame(
  Activity = character(),
  Coverage = numeric(),
  Detection_Other = numeric(),
  stringsAsFactors = FALSE
)

for(act in activities) {
  cat(sprintf("--- Processing Class: %s ---\n", act))
  
  # Force character matching to avoid factor level issues
  idx_act <- which(as.character(window_labels) == act)
  
  if(length(idx_act) < 10) {
    cat(sprintf("Not enough data for %s (Found %d windows). Skipping...\n",
                act, length(idx_act)))
    next
  }
  
  fd_act <- acc_fd[idx_act]
  
  # Split into Train and Calibration
  set.seed(42)
  n_train <- max(5, floor(0.6 * length(idx_act)))
  idx_train <- sample(1:length(idx_act), n_train)
  idx_calib <- setdiff(1:length(idx_act), idx_train)
  
  fd_train <- fd_act[idx_train]
  fd_calib <- fd_act[idx_calib]
  
  # Perform FPCA on Training Data
  nharm_max <- min(10, n_train - 1)
  pca_act <- pca.fd(fd_train, nharm = nharm_max)
  
  # Keep PCs explaining at least 90% variance
  cumvar <- cumsum(pca_act$varprop)
  p_keep <- which(cumvar >= 0.90)[1]
  if (is.na(p_keep)) p_keep <- min(3, nharm_max)
  cat(sprintf("  Using %d PCs (explains %.1f%% var)\n",
              p_keep, cumvar[p_keep]*100))
  
  # FPCA Harmonics Plot
  tryCatch({
    par(mfrow = c(1, min(2, p_keep)))
    plot(pca_act$harmonics[1:min(2, p_keep)],
         xlab = "Time (s)", ylab = "Harmonic")
    mtext(paste("FPCA Perturbations:", act), side = 3,
          line = -1, outer = TRUE, font = 2)
    par(mfrow = c(1, 1))
  }, error = function(e) {
    cat("Warning: Could not plot FPCA harmonics:", e$message, "\n")
  })
  
  # Define conformity scores (Mahalanobis distance)
  scores_train <- pca_act$scores[, 1:p_keep, drop = FALSE]
  mu_scores    <- colMeans(scores_train)
  Sigma_scores <- cov(scores_train)
  
  # Add regularization for numerical stability
  Sigma_scores <- Sigma_scores + diag(1e-6, p_keep)
  
  # Project Calibration data
  scores_calib_raw <- inprod(fd_calib, pca_act$harmonics[1:p_keep])
  if(is.vector(scores_calib_raw)) {
    scores_calib <- matrix(scores_calib_raw, ncol = p_keep, byrow = FALSE)
  } else {
    scores_calib <- scores_calib_raw
  }
  
  # Ensure correct dimensions
  if(ncol(scores_calib) != p_keep) {
    scores_calib <- t(scores_calib)
  }
  
  d2_calib <- mahalanobis(scores_calib, center = mu_scores, cov = Sigma_scores)
  
  # Calculate Threshold (alpha = 0.10)
  q_hat <- as.numeric(quantile(d2_calib, probs = 1 - conf_alpha, type = 1))
  cat(sprintf("  Ellipsoidal Threshold (q_hat): %.3f\n", q_hat))
  
  # Empirical coverage for this activity (calibration windows)
  coverage_act <- mean(d2_calib <= q_hat)
  cat(sprintf("  Empirical coverage for %s (should be ≈ %.2f): %.3f\n",
              act, 1 - conf_alpha, coverage_act))
  
  # how often do other activity windows fall OUTSIDE this band's ellipsoid?
  idx_other_global <- which(window_labels != act)
  if (length(idx_other_global) > 0) {
    fd_other <- acc_fd[idx_other_global]
    
    # Project other windows into this class's PC space
    scores_other_raw <- inprod(fd_other, pca_act$harmonics[1:p_keep])
    if (is.vector(scores_other_raw)) {
      scores_other <- matrix(scores_other_raw, ncol = p_keep, byrow = FALSE)
    } else {
      scores_other <- scores_other_raw
    }
    
    if(ncol(scores_other) != p_keep) {
      scores_other <- t(scores_other)
    }
    
    d2_other <- mahalanobis(scores_other, center = mu_scores,
                            cov = Sigma_scores)
    
    detection_other <- mean(d2_other > q_hat)
    cat(sprintf("  Detection of non-%s windows (d2 > q_hat): %.3f\n",
                act, detection_other))
  } else {
    detection_other <- NA
  }
  
  # Store results
  coverage_results <- rbind(coverage_results,
                            data.frame(Activity = act,
                                       Coverage = coverage_act,
                                       Detection_Other = detection_other))
  
  # Identify outliers for this class in score space
  is_outlier <- d2_calib > q_hat
  cat(sprintf("  Found %d outliers in %s calibration set.\n",
              sum(is_outlier), act))
  
  # Build a data frame with first two PCs for calibration windows
  if (ncol(scores_calib) >= 2) {
    pc1_vals <- scores_calib[, 1]
    pc2_vals <- scores_calib[, 2]
  } else {
    pc1_vals <- scores_calib[, 1]
    pc2_vals <- rep(0, nrow(scores_calib)) # Pad with 0s if only 1 PC
  }
  
  df_scores_calib <- data.frame(
    PC1 = pc1_vals,
    PC2 = pc2_vals,
    Outlier = factor(ifelse(is_outlier, "Outlier", "Inlier"),
                     levels = c("Inlier", "Outlier"))
  )
  
  # Outlier Plot
  p_out <- ggplot(df_scores_calib, aes(x = PC1, y = PC2, color = Outlier)) +
    geom_point(size = 2, alpha = 0.8) +
    scale_color_manual(values = c(Inlier = "grey60", Outlier = "red")) +
    labs(title = sprintf("FPCA Scores for %s (Calibration)", act),
         subtitle = "Outliers defined by Mahalanobis distance > q_hat",
         x = "PC1 score", y = "PC2 score") +
    theme_minimal() +
    theme(legend.position = "top",
          panel.grid.major = element_line(color = "grey90"))
  
  print(p_out)
  
  # Build Prediction Band for Time Domain Plotting
  phi_mat <- eval.fd(grid_eval, pca_act$harmonics[1:p_keep])
  if(is.vector(phi_mat)) phi_mat <- matrix(phi_mat, ncol = p_keep)
  
  mean_vec  <- as.vector(eval.fd(grid_eval, pca_act$meanfd))
  
  band_center <- numeric(length(grid_eval))
  band_upper  <- numeric(length(grid_eval))
  band_lower  <- numeric(length(grid_eval))
  
  for (r in seq_along(grid_eval)) {
    phi_vec <- phi_mat[r, ]
    if(length(phi_vec) != p_keep) phi_vec <- matrix(phi_vec, nrow=1)
    
    center_t <- mean_vec[r] + sum(mu_scores * phi_vec)
    var_t <- as.numeric(phi_vec %*% Sigma_scores %*% t(phi_vec))
    radius_t <- sqrt(max(q_hat * var_t, 0))
    
    band_center[r] <- center_t
    band_upper[r]  <- center_t + radius_t
    band_lower[r]  <- center_t - radius_t
  }
  
  # Store band data
  df_bands_all <- rbind(df_bands_all, data.frame(Activity = act,
                                                 t = grid_eval,
                                                 center = band_center,
                                                 upper = band_upper,
                                                 lower = band_lower))
  
  # Reconstruct calibration curves for plotting within the band
  proj_calib <- matrix(NA, nrow = length(grid_eval), ncol = nrow(scores_calib))
  for (i in 1:nrow(scores_calib)) {
    proj_calib[, i] <- mean_vec + as.numeric(phi_mat %*% scores_calib[i, ])
  }
  
  df_proj_all <- rbind(df_proj_all, data.frame(
    Activity = act,
    t = rep(grid_eval, ncol(proj_calib)),
    val = as.vector(proj_calib),
    curve = rep(paste0(act, "_", 1:ncol(proj_calib)), each = length(grid_eval))
  ))
}
## --- Processing Class: Sedentary ---
##   Using 1 PCs (explains 99.8% var)

##   Ellipsoidal Threshold (q_hat): 1.751
##   Empirical coverage for Sedentary (should be ≈ 0.90): 0.901
##   Detection of non-Sedentary windows (d2 > q_hat): 0.933
##   Found 16 outliers in Sedentary calibration set.

## --- Processing Class: Walking ---
##   Using 1 PCs (explains 99.9% var)

##   Ellipsoidal Threshold (q_hat): 5.841
##   Empirical coverage for Walking (should be ≈ 0.90): 0.906
##   Detection of non-Walking windows (d2 > q_hat): 0.041
##   Found 8 outliers in Walking calibration set.

## --- Processing Class: Running ---
##   Using 1 PCs (explains 90.9% var)

##   Ellipsoidal Threshold (q_hat): 6.175
##   Empirical coverage for Running (should be ≈ 0.90): 0.900
##   Detection of non-Running windows (d2 > q_hat): 0.060
##   Found 18 outliers in Running calibration set.

## Coverage results:
##    Activity  Coverage Detection_Other
## 1 Sedentary 0.9012346      0.93323217
## 2   Walking 0.9058824      0.04103165
## 3   Running 0.9000000      0.06006494

The coverage table shows that all three conformal sets achieve empirical coverage essentially equal to the nominal 90%, confirming that the inductive conformal construction in FPCA score space works as intended for each activity. Sedentary also has a very high Detection_Other (~0.93), meaning almost all non sedentary windows fall outside the sedentary ellipsoid, which fits the strong separation seen in the mean/boxplot analyses. Walking and Running, instead, have low Detection_Other (~0.04–0.06), indicating that many non walking and non running windows still lie inside their ellipsoids; this is consistent with their higher variability and partial overlap in latent space.

The class specific FPCA score plots make this explicit: outliers for Running and Walking sit on the right tail of PC1 but still form a continuum with inliers, while Sedentary outliers are almost all separated in PC1.

Finally, the single retained harmonic per class explains >90–99% of variance and has smooth, low frequency shape, confirming that a one dimensional Gaussian model in score space is a good approximation and that ellipsoidal conformal sets are a natural choice for these activity specific functional signatures.

Conformal Interval Visualization

Here we visualize the boundaries of “normal” behavior for each activity and the actual visual proof of our conformal prediction model.

To illustrate these boundaries, our first key visualization is the Time-Domain Conformal Bands plot . This translates the abstract ellipsoidal math back into the physical world. Specifically, it plots the 90% confidence “tube” (ribbon) over time for each activity. The practical application here is straightforward: if a new 5 second window of movement is recorded, and its curve fits completely inside the “Walking” tube, the model confidently classifies it as walking.

While the first plot focuses on the temporal aspect, the Global Latent Space Ellipsoids plot provides a broader, structural perspective. Instead of looking at curves over time, this plot zooms out and examines the overall mathematical separation of the activities. To achieve this, it runs a global Functional PCA to compress all curves into just two principal variables (PC1 and PC2). Finally, it plots every single 5 second window as a dot, drawing 90% confidence boundaries (ellipses) around the clusters to clearly demonstrate how well the Sedentary, Walking, and Running states separate from each other in the latent space.

if (nrow(df_bands_all) == 0) {
  cat("\n[ERROR] No data was generated for the bands\n")
} else {
  
  # Ensure factors for consistent plotting order
  df_bands_all$Activity <- factor(df_bands_all$Activity, levels = activities)
  df_proj_all$Activity <- factor(df_proj_all$Activity, levels = activities)
  
  # Conformal Bands
  p_bands_multi <- ggplot() +
    geom_ribbon(data = df_bands_all, aes(x = t, ymin = lower, ymax = upper, fill = Activity),
                alpha = 0.3) +
    geom_line(data = df_proj_all, aes(x = t, y = val, group = curve, color = Activity),
              alpha = 0.4, linewidth = 0.4) +
    geom_line(data = df_bands_all, aes(x = t, y = center, color = Activity),
              linewidth = 0.8, linetype = "dashed") +
    facet_wrap(~Activity, scales = "free_y") +
    scale_fill_manual(values = col_map) +
    scale_color_manual(values = col_map) +
    labs(title = "Conformal Prediction Bands for All Activities (90% Coverage)",
         subtitle = "Gaussian Mixture Components: Projections from FPCA Mahalanobis Ellipsoids",
         x = "Relative Time (seconds)", y = "Acceleration Magnitude") +
    theme_bw() +
    theme(legend.position = "none",
          strip.background = element_rect(fill="grey90"),
          strip.text = element_text(face="bold", size=12),
          panel.grid.major = element_line(color = "grey95"))
  
  print(p_bands_multi)}

# Global Latent Space Conformal Ellipsoids
pca_global <- pca.fd(acc_fd, nharm = 2)
df_latent <- data.frame(
    PC1 = pca_global$scores[, 1],
    PC2 = pca_global$scores[, 2],
    Activity = factor(window_labels, levels = c("Sedentary", "Walking", "Running")))
df_latent <- na.omit(df_latent)

get_ellipse_points <- function(data, level = 0.90, npoints = 100) {
  center <- colMeans(data[, c("PC1", "PC2")])
  cov_mat <- cov(data[, c("PC1", "PC2")])
  chi_val <- qchisq(level, df = 2)
  eig <- eigen(cov_mat)
  angles <- seq(0, 2*pi, length.out = npoints)
  ellipse <- t(center + sqrt(chi_val) * eig$vectors %*% 
               diag(sqrt(eig$values)) %*% 
               rbind(cos(angles), sin(angles)))
  
  data.frame(x = ellipse[, 1], y = ellipse[, 2])
}

p <- plot_ly() %>%
  layout(title = "Latent Space Conformal Regions (90% Coverage)",
         xaxis = list(title = "PC1"),
         yaxis = list(title = "PC2"),
         dragmode = "zoom")

for(act in unique(df_latent$Activity)) {
  df_act <- df_latent[df_latent$Activity == act, ]
  
  if(nrow(df_act) > 5) {
    ell_pts <- get_ellipse_points(df_act, level = 0.90)
    
    
    p <- p %>%
      add_trace(x = ell_pts$x, y = ell_pts$y,
                type = 'scatter', mode = 'none',
                fill = 'toself',
                fillcolor = col_map[act],
                opacity = 0.15,
                line = list(width = 0),
                name = paste(act, "region"),
                showlegend = FALSE,
                hoverinfo = 'none')
    
    
    line_width <- ifelse(act == "Sedentary", 10, 3)
    
    p <- p %>%
      add_trace(x = ell_pts$x, y = ell_pts$y,
                type = 'scatter', mode = 'lines',
                line = list(
                  color = col_map[act], 
                  width = line_width,
                  shape = "spline"
                ),
                name = paste(act, "ellipsoid"),
                showlegend = TRUE,
                hoverinfo = 'none')
  }
  
  
  point_size <- ifelse(act == "Sedentary", 7, 5)
  
  p <- p %>%
    add_trace(data = df_act,
              x = ~PC1, y = ~PC2,
              type = 'scatter', mode = 'markers',
              marker = list(
                size = point_size, 
                opacity = 0.8, 
                color = col_map[act]),
              name = act,
              hoverinfo = 'text',
              text = ~paste(act, "<br>PC1:", round(PC1, 2), "<br>PC2:", round(PC2, 2)))}

p

The latent space plot shows all 5 second windows projected onto the first two FPCA components, with separate Gaussian ellipses at 90% level for Sedentary, Walking, and Running. Sedentary points form a short, tight segment near the origin, while Running spreads widely along PC1 and PC2, illustrating why Sedentary is easiest to separate in conformal score space. Walking sits between them, overlapping partially with both, which explains its lower Detection_Other in the coverage table.

The time domain conformal bands panel then shows each class ellipsoid into acceleration space with narrow, low amplitude tubes for Sedentary, and progressively higher and wider tubes for Walking and Running. Most projected calibration curves stay inside their respective tubes, visually confirming the ~90% empirical coverage and showing that the FPCA Mahalanobis conformal approach captures both intensity and variability of each activity.

Conclusions

The application of Conformal Prediction based on Functional PCA (FPCA) and Mahalanobis distance provided highly robust results for modeling different movement intensities in our accelerometer data. Extracting just a single Principal Component was sufficient to explain nearly all the variance for each class: 99.8% for Sedentary, 99.9% for Walking, and 90.9% for Running. Building on this robust dimensionality reduction, when setting a target confidence level of 90%, the algorithm’s empirical coverage on the calibration and test data hits the nominal target almost perfectly. We recorded coverages of 90.1% for Sedentary, 90.6% for Walking, and 90.0% for Running. This mathematically demonstrates that the conformal regions, built by calculating the threshold on the Mahalanobis distances, are perfectly calibrated. As a result, the conformal tubes in the time domain enclose exactly the desired proportion of trajectories for each class, ensuring the distribution free validity of the method.

Furthermore, analyzing the curves classified as anomalies alongside the Detection_Other metric, which measures the percentage of windows from other classes correctly rejected by a given class’s ellipsoid, reveals a fundamental geometric property of human movement within the latent space. The resting, Sedentary class generates an extremely tight ellipsoid. This allows the model to perfectly isolate micro movements and reject almost everything that is not rest, resulting in an exceptionally high Detection_Other rate of 93.3%. The 16 outliers found within this class during calibration likely represent slight bumps, posture changes, or sensor jitter. Conversely, Walking and Running activities possess significantly higher energy and variance, translating into much larger ellipsoids with values of 5.841 and 6.175, respectively. From a geometric perspective, these high variance ellipsoids inevitably encompass or “swallow” very low variance curves. Consequently, the conformal bands for Walking and Running naturally contain the Sedentary curves, leading to much lower true negative rejection rates, which is reflected in a Detection_Other of 4.1% for Walking and 6.0% for Running.